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SUMMARY 


The Earth Radiation Budget Experiment (ERBE) has given rise to numerous analyses 
involving various data sets. Spatial coordinate transformations are developed for 
several different grid systems used for the compressed storage of global satellite 
data. These include grids for the ERBE project, for analysis of existing satellite 
radiation data, and for the U.S. Air Force cloud data base. Generally the grids are 
defined in terms of a two-dimensional index system, with associated regional box num- 
bers, The transformations associate longitude and latitude with a particular grid 
location, and vice versa. In addition, the relationship between longitude and right 
ascension is given. This transformation combines position information for the Sun 
with longitude-latitude coordinates as the Earth rotates in inertial space; from it, 
certain information about the Sun's position relative to an observer can be obtained. 
The transformations developed in this memorandum are implemented in FORTRAN IV in an 
appendix. Such implementation is suggested by the fact that the spatial grid systems 
are defined, rather than just described, by their transformation algorithms and that 
actual computer- coded versions of the transformations are the best way to avoid some 
potential ambiguities in each system. 


INTRODUCTION 

Development of the Earth Radiation Budget Experiment (ERBE, described in ref. 1) 
has resulted in a large number of directly related and supporting software tasks. 

The latter involve other existing or potential data sources, including archived U.S. 
Air Force cloud data (the so-called 3DNEPH data) , radiation data from the eastern 
Geostationary Operational Environmental Satellite (GOES-east, hereinafter referred to 
simply as GOES) , and the Nimbus-Earth Radiation Budget experiment. Examination of 
such data requires many spatial transformations back and forth between locations on 
the Earth's surface (longitude and latitude) and a grid system suitable for analysis 
and compressed storage of large quantities of information. In addition, analysis of 
variables which exhibit diurnal behavior often requires temporal transformations 
between a standardized global time and a local clock or solar time. Thus, this memo- 
randum has two purposes: first, to provide spatial transformation equations relating 

position on the Earth's surface to locations in a grid system and, second, to provide 
transformations relating inertial space to coordinate systems on the Earth's rotating 
surface. The latter equations are equivalent to relating universal time (UT) to a 
local or solar time; they allow position information for the Sun and an observer to 
be compared in the same coordinate system. 

Spatial transformations in algorithm form not only make data manipulation pos- 
sible but also serve to actually define the coordinate systems involved. From a 
computational point of view, careful definitions are especially important at the 
boundaries of grid boxes. One way of looking at this problem is to ask whether points 
exactly on the Equator belong in Northern or Southern Hemisphere boxes. Similar care 
is needed at each pole and at a longitude of 0° and 360°, where indexing algorithms 
may produce indices outside the expected range. Pictures of the grid do not resolve 
these details, so that precise algorithms are necessary to avoid ambiguity and compu- 
tational problems in a computer code. 



SYMBOLS 


a relative azimuth between an observer and the Sun, as defined in figure 6, 

deg 

A array used in defining the Nimbus-ERB grid 

B box number in any of the grid systems 

C variable related to longitude, for defining the Nimbus-ERB grid, deg 

COL matrix colimrn identifier within boxes of the 3DNEPH grid 

d integer part of ~ l)/8, for defining the 3DNEPH grid 

D array used in defining the GOES grid 

f angular variable used in defining the 3DNEPH grid 

G array used in defining the Nimbus-ERB grid 

I longitude increment in a latitude band of the Nimbus-ERB grid, deg 

INT FORTRAN function: INT (m) = (sign m) (Largest integer ^ |m|) 

L geocentric longitude, measured positively to the east from the Greenwich 

meridian, deg 

LAT integer latitude band identifier 

LON integer longitude box identifier 

m dummy variable for defining INT and MOD functions 

MOD FORTRAN function: MOD(m,n) = m - n INT(m/n) 

n number of boxes in a latitude band of an ERBE grid, or dummy variable for 

defining the MOD function 

r position vector, any units 

ROW matrix row identifier within boxes of the 3DNEPH grid 

s size of sides on one of three ERBE grids, 2.5°, 5°, or 10° (also see 

subscript list) 

T time in Julian centuries (36 525 days) from Jan. 0.5, 1900, Julian Date 

241 5020.0 

UT universal time, hours, minutes, and seconds, or fractions of a day 

u,v dimensionless variables used in defining the 3DNEPH grid 

x,y,z vector components, any units 
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Z array used in defining the Langley GOES grid 

a right ascension, deg 

6 declination, deg 

£ angle between ecliptic plane and Earth's Equator, =23.45° 

0 zenith angle, deg 

\ geocentric latitude, deg (90° - X - -90°) 

geocentric colatitude, deg (0° ^ A^^ ^ 180°) 

T first point of Aries 

to Earth's rotation rate, 360.9856473 deg/day 

Subscripts : 

E Earth 

g Greenwich meridian 

g, o Greenwich meridian at o'^ UT 

grd a place on the Earth's surface 

G Langley GOES grid 

m meridian 

mid midpoint of grid box 

min minimum 

max maximum 

N Nimbus-ERB grid 

obs orbiting observer of the Earth's surface 

s one of the three ERBE grid systems when used with LAT, LON, or B; or to 

the Sun in the discussion of the relationship between longitude and right 
ascension 

3D 3DNEPH grid 

Superscripts : 

h, m,s hours, minutes, seconds 

dimensioned vector 
unit vector 
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THE ERBE GRID SYSTEM 


Selection of a suitable grid system for the ERBE project is an important process, 
affecting literally millions of computer calculations during the life of the project. 
Representing spacecraft data on a grid system implies that it makes sense to portray 
information at the spatial resolution of the grid. A basic requirement of the ERBE 
project is that it produce "regional" monthly averages of the Earth's radiation budget. 
A "region" is loosely defined as an area a few hundred (250-500) kilometers on a side 
(ref. 1) . Within the constraints of this requirement, several different grid systems 
have been proposed. "Constant area" systems start with a box of appropriate size at 
the Equator and seek to define longitude-latitude boundaries over the globe in a 
systematic way which preserves at least approximately the area of the equatorial ref- 
erence box. This is done in one of two ways: either the latitude boundaries of the 

box are kept constant while the longitude boundaries are adjusted, or vice versa. An 
alternative system gives up constant area in favor of constant longitude-latitude 
boundaries. This leads to boxes which vary greatly in area from Equator to pole. 


The ERBE Science and Data Management Teams have informally evaluated several pro- 
posed grid systems on the basis of meeting resolution requirements, ease of nesting 
from one resolution to another, ease of indexing, and possibilities for a reasonable 
global representation of the ERBE data products. The equal longitude-latitude system 
has been chosen over any of the equal area systems, despite the problems of inter- 
preting regional data in grid boxes which vary in area by more than an order of magni- 
tude. The layout of a 2.5° system is given in figure 1. To simplify indexing, lati- 
tude A is replaced with colatitude, where A^q = 90° - A, so that 0 ° ^ ^ 180°. 

The boxes at the Equator have areas of 77 426 km^, while those nearest the poles have 
areas of 1686 km^. Because of its simple boundary definition (fixed increments in 
longitude and colatitude) , the ERBE grid is easy to define and work with. Actually, 
three different grid systems are going to be employed in ERBE analysis, with angular 
resolutions s of 2.5°, 5°, and 10°. These contain 10 368, 2592, and 648 boxes, 
respectively, the number of boxes n in any latitude band being 144, 72, and 36, 
respectively. The following equations define each of the systems by giving trans- 
formations from spatial coordinates to a location in the grid, and vice versa. 
Transformations between the systems are also given; all the transformations in this 
memorandum are defined specifically in terms of geocentric, as opposed to geodetic, 
coordinates. Given a colatitude A^q and longitude L (0° ^ L ^ 360°) , the location 
in the grid (LAT^,LONg) is 


LAT3 


|lNT(A^ys + 1) 


n /2 


(Aco ^ 180 °) 
(Aco = 180°) 


LON 

s 


INT(L/s + 1) 


(L 7 ^ 360°) 
(L = 360°) 


The INT function truncates (does not round off) its argument: 


INT(m) = (sign m) ^Largest integer - |m|) 
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The specific definitions at = 180° and L = 360° are given to prevent, for 

example, LAT 2 5 = 73 and LON 2 5 = 145. Each test for these limits in computer- 
coded form requires the equivalent of a FORTRAN IF statement and so should be avoided 
if possible by restraining the values of X^^ and L. This grid system numbers the 
latitude boxes from north to south and longitude boxes from west to east. 

A consecutively numbered box system starting at the North Pole and Greenwich 
meridian, and winding eastward around the globe is given by 


Bg = (LATg - l)n + LONg 


Note that this system is defined so that in the 2.5° grid, for example, a point at 
X^^ = 90° and L = 0° is in the upper left-hand corner of the box identified by 
LATg = 37 and LONg = 1, or Bg = 5185. That is, points on the Equator belong in 
Southern Hemisphere boxes. More generally, each box contains points on its northern 
and western sides, but not on its southern and eastern sides. Given a box Bg, LATg, 
and LONg are given by 


LATg = INT[ (Bg - l)/n] + 1 
LONs = Bg - (LATg - l)n 


Matrix references at one resolution can be obtained from those at another. Let 
Xg = LATg or LONg. Then 

Xg = INT[ (X 2 5 - l)/2] + 1 
Xj_Q = INT[ (X 2,5 - D/4] + 1 
= INT[ (X 5 - l)/2] + 1 


Consecutively numbered boxes at different resolutions can be related as follows: 


Bg = 72 INT [(B 2^5 - D/288] + INt{mOD [( 62^5 - l),144]/2} + 1 
Bj^q = 36 INT [(B 2_5 - D/576] + INt{mOD [(B 2 5 - l),144]/4^ + 1 
= 36 INT [(Bg - D/144] + INt{mOD [(Bg - l),72]/2} + 1 


v/here MOD(m,n) = m - n INT(m/n). The reverse calculation finds all the smaller 
boxes included in a given larger box. First find the upper left-hand corner box: 

B 2 . 5 (corner) = 576 INt[ (B^^q - l)/36] + 4 M0D[ (B^^q - 1),36] + 1 
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The remaining 15 boxes in the larger box are 


B 2 _ 5 (corner) + 1,2,3,144,145,146,147,288,289,290,291,432,433,434,435 

Similarly, 

85 (corner) = 144 INT [(B^q - l)/36] + 2 MOD [(B^^q - 1),36] + 1 
B^(corner) + 1,72,73 

and 

B 2 , 5 (corner) = 288 INT[ (B 5 - l)/72] + 2 M0D[ (B 5 - 1),72] + 1 
B 2 ^(corner) + 1,144,145 

Given a particular box identified by LONg and LATg, the longitude and colatitude 
coordinates of the midpoint are given by 

^co,mid “ “ l)s + s /2 

^mid = - l)s + s /2 


OTHER GRID SYSTEMS OF INTEREST 


The Langley GOES Grid 

ERBE-related research at NASA Langley Research Center has included significant 
use of GOES data. This has been processed at a basic resolution of 8 km and finally 
on a constant area grid system which starts with 2.25° x 2.25° boxes at the Equator 
(ref. 2) . Much of the work has involved 1600 grid boxes in the Western Hemisphere, 
centered around the Equator and extending to ±45° in latitude. The Langley GOES grid 
is shown in figure 2. The conversion from latitude and longitude to a 2-digit refer- 
ence in the GOES grid system is as follows: 


LATg 


[iNT[ (47.25 - A)/2.25] 
I 40 


(-45° < A < 45°) 
(A = -45°) 


LON^ = INT[(Z(i) + L)/D(i)] 
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where 


i = INt(|X|/ 18 + l) 

and 


i 

Z(i) 

D(i) 

^in 


1 

123.75 

2.25 

-121.5 

-31.5 

2 

132.50 

2.50 

-130.0 

o 

o 

3 

141.00 

3.00 

-138.0 

0 

00 

1 — ( 

1 


- L^in- 


The GOES grid boxes are numbered consecutively from top (+45° latitude) to bottom 
(-45°) , west to east, with 40 boxes in each row. Each box contains points on its 
northern and western sides, but not on its southern and eastern sides, so that points 
directly on the Equator belong in Southern Hemisphere boxes. The box number is 

Bg = (LATj^ - 1)40 + LONq 

The reverse procedure is to find the longitude and latitude of the midpoint of each 
GOES box. Given a box number Bq, 

LATq = INT[ (Bq - l)/40] + 1 

LONq = Bq - (LATq - 1)40 

The midpoints are 


Vd = 46.125° - 2.25 °LATq 

Lmid = D(i) LONq - [Z(i) - D(i)/2] 

where 


i = INT[ (28 - k)/8] 


k = 


[lATq 

(IATq < 

20) 

^1 - LATq 

(IATq > 

20) 
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The Nimbus-ERB Grid System 


The Nimbus Earth Radiation Budget (ERB) Experiment has been conducted on the 
Nimbus 6 and 7 spacecraft (ref. 3) . Radiation budget data have been analyzed globally 
in "target areas" defined by a constant area grid system of the type considered, but 
not chosen, for the ERBE. It consists of 2070 consecutively numbered boxes starting 
at the South Pole and spiralling westward up to the North Pole. Each box contains 
points on its southern and eastern sides, but not on its northern and western sides, 
so that points directly on the Equator belong in Northern Hemisphere boxes. (Note 
that this is opposite from the ERBE and Langley GOES grids.) The grid system is 
pictured in figure 3. Boxes at the Equator are 4.5° x 4.5°; they are defined as 
follows : ^ 


LATj^ = 


'iNT[ (A + 90°) /4. 5] + 1 
140 


(-90° < A < 90°) 
(A = 90°) 


LONj^ = INT(C/I) + 1 


where for 0° - L - 360°, 
C = 360° - L 


and the longitude increment for each value of LATj^ is 


I = 360°/A(i) 



(LATj^ < 20) 
(LATjj > 20) 


The array A gives the number of boxes in each latitude band from pole to Equator. 
(The same niomber in A applies to corresponding latitude bands in each hemisphere.) 
Then 


G(i) + LON 




N 


[G(i) + LONj^ - A(i) 


(L 0) 
(L = 0) 


The array G gives the sequential number of the last box in the previous latitude 
band, starting at the South Pole. The arrays A and G are 


^These equations are based on FORTRAN subroutines supplied by Guy T. Germana, 
Systems and Applied Sciences Corporation, Riverdale, Maryland, March 1981. 
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A: 

3,9,16 

,20,30 

,36,40 

,45,48 

,60,60 

,60,72 

,72,72 

,72, 

G: 

0 

3 

12 

28 

48 

78 

114 

154 


199 

247 

307 

367 

427 

499 

571 

643 


715 

795 

875 

955 

1035 

1115 

1195 

1275 


1355 

1427 

1499 

1571 

1643 

1703 

1763 

1823 


1871 

1916 

1956 

1992 

2022 

2042 

2058 

2067 


The reverse process gives midpoint longitude and latitude for each box. First 
find the smallest element j in array G such that G(j) - Then 


LAT 


N 



1 


(1 5 Bn ^ 2067) 
(2068 ^ Bn ^ 2070) 


The column identifier is 


LON^ = [b^ - G(LATj^)] 


For I, as previously defined, the midpoint longitude and latitude are 


= 360° - [KLONn - 1) + 1/2] 
/mid = "90° + (LATj^ - 1)4.5° + 2.25° 


The 3DNEPH Grid System 


For a number of years the U.S. Air Force has produced global cloud cover infor- 
mation which is archived at the Environmental Technical Applications Center, Ashe- 
ville, North Carolina. This is available to civilian users through the National 
Climatic Center, Asheville, North Carolina. The grid system is an 8 x 8 set of boxes 
superimposed on a polar stereographic projection, one set of 64 boxes for each 
hemisphere. Points on the Equator belong in Northern Hemisphere boxes. This grid is 
described in reference 4 and pictured in figure 4. The four corner boxes are not 
used. Alignment of the grid along L = 10° is an apparently arbitrary choice. Each 
box is subdivided into a 64 x 64 array of points spaced a nominal 25 nautical miles 
apart. Point (1,1) occupies the upper left-hand corner of its box as shown in fig- 
ure 4(c). Longitude (-180° ^ L ^ 180°) and latitude are related to a particular box 
and to a row and column in that box as follows: 


u 


cos (350° + L) 
cos (350° + L) 


(62317.6272) (l - sinlX|) 


(l + sin|X|)[l + tan^(350° + L) 


(L / 100° or -80°) 
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Since u is undefined for L = 100° or -80°, it is necessary (and sufficient) for 
computation to let 


L = 99.9999° 

( 1 100° - L 1 

< 10 

L = -79.9999° 

(1-80° - l| 

< 10“'^) 


Then 


V = u tan (350° + L) 


Let 


i = 257 + INT(u + 0.5) 
j = 257 - sign (A) [ INT (v + 0.5)] 


Then let 


ROW = j - 64 INT [(k - l)/8] 

where k is the smallest nonnegative integer such that 1 ^ ROW ^ 64 
COL = i - 64£ 


where £ is the smallest nonnegative integer such that 1 - COL - 64 
and COL are the row and column identifiers in the box numbered 


As an example of locations on the 3DNEPH grid, the grid point closest 
and 0° latitude is in Northern Hemisphere box 40, row 43, column 55. 
operation is to convert ®3D' ROW, and COL into longitude and latitude 


L = MOD [(tan“l (v/u) - 350°), 360°] 


A 


sin~^f 


-sin“^f 


(Northern Hemisphere) 
(Southern Hemisphere) 


where 

sin f = [62317.6272 - (u^ + v2) ]/[62317.6272 + (u2 + v2) ] 


u = i - 257 


Further, let 


Finally, ROW 


to 0° longitude 
The reverse 
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V 


257 - j 


(Northern Hemisphere) 


j - 257 


(Southern Hemisphere) 


i = [ (B 3 P - 1) - 8d]64 + COL 
j = 64d + ROW 
d = INT [( 63 ^^ - l)/ 8 ] 


Care must be taken to obtain the proper quadrant from the inverse tangent function. 
The longitude and latitude of the grid point for the above example (box 40, row 43, 
column 55) are L = 0.311° and X = 0.017°. 


THE RELATIONSHIP BETWEEN LONGITUDE AND RIGHT ASCENSION 

The astronomical equations relating the Sun's position to the Earth are given in 
a (nearly) inertial framework; an Earth-equatorial system with the x-axis pointing 
to a particular place in the sky (the first point of Aries, denoted by T) which 
remains fixed for most practical purposes. The basic time system for processing ERBE 
data is universal time (UT) . This is a globally constant mean solar time in which 

O^o'^O^ UT starts the calendar day (refs. 5 and 6 ). UT has the advantage of being a 
centralized time system which does not depend on an observer's location on the globe. 
However, there are disadvantages when data must be presented which show diurnal varia- 
tions related to solar activity. If 12*^ UT is local noon at one place on the globe, 
it is midnight on the other side of the globe. Physically, then, analyses of diurnal 
behavior are better treated in a local time system with, perhaps, a UT time scale also 
maintained to retain a standardized link to other data. The equivalence between local 
time and UT is in one sense a spatial transformation between the inertial coordinate 
system and the Earth's own rotating longitude-latitude system. UT is defined from the 
Greenwich meridian (0° longitude), so that it is local mean solar time ("clock time") 
at Greenwich ± 7.5° of longitude. Every 15° of longitude corresponds to a 1-hour 
shift, neglecting local variations and "daylight savings" time. The adjustment is 
measured positively to the east, so that moving east by 15° requires setting a clock 
ahead by 1 hour. Given a particular UT and measuring longitude positively from 
Greenwich eastward. 


Local hour 


MOD 


(^UT hour + INT [(L + 7.5°)/15°]^ 



Whereas longitude is the measure of angular distance in the Earth's rotating 
coordinate system, right ascension is the corresponding measure in the inertial 
system. Longitude and right ascension are related through the equation of time, which 
contains astronomical information about the Earth's movement relative to the Sun. 

The geometry is illustrated in figure 5. At Greenwich midnight, the start of a 
calendar day, the right ascension of the Greenwich meridian is 
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^ = 99.6909833° + 36000. 7689°T + 0 . 00038708°T^ 
y / u 

where T is the time in Julian centuries from January 0.5, 1900: 


T = (Julian date - 241 5020.0)/36525 


A table of Julian dates and a detailed explanation of the use of these equations can 
be found in reference 6. To determine the right ascension (called sidereal time) of 
any meridian on the Earth's surface. 


a 


m 


cx_ ^ + (jJt + L 
g ,o m 


where t is fraction of a day measured from O^O^O^ UT, co is the Earth's rotation 
rate, 360.9856473 deg/day, and is the longitude of the meridian. 

THE RELATIONSHIP BETWEEN AN OBSERVER AND THE SUN 

The Sun's location in geocentric inertial space can be determined from the 
equations of the Earth's orbital motion around the Sun. Let rg be the Earth's 
position vector in heliocentric ecliptic space, where the ecliptic plane is defined 
as the plane in which the Earth travels around the Sun. (The necessary equations can 
be found in reference 5 or 6, or in the implementation of this transformation in the 
appendix.) Then, the unit vector defining the Sun's position relative to the Earth 
in this same coordinate system is 


^S = (-^e/I^eI'“^e/I^eI'"^e/I^eI) 


The Sun's unit position vector is rotated from the ecliptic to an Earth-equatorial 
system: 

^s = ^s 

Ys = Ys cos e - Zg sin £ 

Zg = Ys sin e + Zg cos e 


where e is the angle between the ecliptic plane and the plane of the Earth's 
equator, about 23.45°. Now, observers on the Earth can relate themselves to the Sun 
by comparing their right ascension and declination to that of the Sun. Declination 
is equivalent to latitude and right ascension is found from the equations in the 
previous section. The geometry is shown in figure 6. In terms of solar right 
ascension Ug and declination 6g, the unit vector of the Sun is 
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Xg = cos 6 ^ cos 

Yg = cos 6 g sin 

Zs = sin 63 

Then 

63 = sin"^Zg 
= tan“^ 


where care must be taken to obtain the proper quadrant from the inverse tangent func- 
tion. The angular distance between the observer's and the Sun's right ascensions is 
a measure of true local time. Noon is defined as the instant at which the declina- 
tions are the same, i.e., when the two meridians coincide. This time is not, in 
general, the same as clock noon because clock time is based on a fictitious mean Sun 
which moves with constant speed along the Equator's projection on the celestial 
sphere. True noon and clock noon always differ by less than 20 minutes, neglecting 
irregularities in local time zones. 

Other useful relationships between an observer and the Sun follow when both 
position vectors are expressed in the same coordinate system. The solar zenith angle 
at a ground observation point, as defined in figure 6 , is 


0s = cos ^(?grd • ^s) 

✓s /s 

If an orbiting observer at ^obs views a point at ^g^d' relative azimuth 

between Sun and observer is defined in figure 7 as 

a = cos-l[(?g^^ ^ ^s) • (Vd ^ ^obs)] 

CONCLUDING REMARKS 

Spatial transformation algorithms have been given for each of four large data 
systems: the Earth Radiation Budget Experiment (ERBE) , eastern Geostationary 

Operational Environmental Satellite (GOES) , Nimbus-Earth Radiation Budget experiment, 
and U.S. Air Force cloud data (called 3DNEPH) . In each case the algorithms deal only 
with data at a particular resolution level, from tens to hundreds of kilometers, which 
is much coarser than the basic resolution of the various systems. The information is 
almost entirely satellite derived from sensors having resolutions of a few kilometers 
(the 3DNEPH cloud data contain other types of observations as well) and has been 
extensively preprocessed before this stage of interest. Thus, the equations presented 
here do not deal directly with the satellite systems themselves, but rather with spe- 
cific, well-documented data products of interest to ERBE scientific investigations. 
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Intercomparison of data between grid systems is not a straightforward matter, and 
this memorandum does not address the large-scale transfer of information from one 
grid to another. Single points may be transformed from one grid to longitude- 
latitude and then onto another grid, but this is not the same thing as expressing the 
distribution of a variable in one system in terms of another system. The transforma- 
tion between longitude and right ascension uses standard astronomical equations; it 
is included for a certain completeness of use to the ERBE project. In particular, 
diurnal models of radiative behavior require specification of solar zenith angle, 
which always necessitates definition of Earth coordinates relative to the Sun. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
November 24, 1981 
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APPENDIX 


FORTRAN IV LISTINGS 

The transformations included in the preceding sections have been implemented in 
FORTRAN IV. A listing of three main programs with subroutines is given below; they 
operate in an interactive mode with NASA Langley Research Center central computers, 
through a remote terminal. The grid system program is LLGRID. There are two sub- 
routines for each grid; one for converting longitude and latitude to grid coordinates 
and another for converting back again. LLGRID expects terminal keyboard input in free 
format (list directed READ) of longitude and latitude, from which it obtains grid 
coordinates for each system. The reverse transformation gives the midpoint longitude 
and latitude of the grid box or, in the case of 3DNEPH, the longitude and latitude of 
the closest grid point. Program LONGRA receives a UT calendar date input from the 
keyboard and converts it to Julian date. It then calculates the right ascension of 
the input longitude. The relationship between a ground-based observer and the Sun is 
illustrated by program ZENITH, which calculates the solar zenith angle for a given 
longitude, latitude, and UT calendar date. Sample input and output follow each 
program listing. 
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APPENDIX 


PROGRAM LLGRID( IN PUT» OUTPUT j TAPES* INPUT# TAPE6* OUTPUT) 
CONVERTS LON AND LAT TO INDICES ON A 3DNE PH#GOE S# NI MBUS# 

OR ERBE GRID. THEN CONVERTS BACK TO CLOSEST GRID LONGITUDE 
AND LATITUDE FOR 3DNEPH GRID AND TO MIDPOINT LONGITUDE AND 
LATITUDE FOR THE OTHER GRIDS. 

DIMENSION IROW(3)# IC0L(3)#IB0X(3) 

C INPUT LONGITUDE AND LATITUDE 

1 PRINTlOlO 

READ*#XLON#XLAT 
IF(XL0N.EQ.999. ) G0T09999 
PRINT1017#XL0N,XLAT 
C 3DHEPH SYSTEM 

CALL LLNEPH{XLQN#XLAT#NEPHB#NROM#NCOL) 

PR INTIOOO 

IF <XLAT.GE.O. ) PRINT1005 
IH»1 

IF(XLAT.LT.O. )IH-2 
IF (XLAT.LT.O. ) PRINT1006 
PRINT1001#NEPHB»NR0W#NC0L 
C ERBE GRID SYSTEM 

PRINT1015 
S»1.2S 
DO 10 I-l>3 
S»S+2. 

CALL LLERB(XLON#XLAT»S# IBOX(I)#IRQW(I)#ICOL( I) ) 

PRINT1003,S 

PRINT1001#IB0X( I ), IROWCI )#IC0UI ) 

10 CONTINUE 

C LONGITUDE# LATITUDE TO GOES GRID 

CALL LLGOES(XLON#XLAT# IGQES#LATG#LONG#IERR) 

IF(IERR.EQ.O) G0T018 

PRINT1007 

G0T019 

18 PRINT1008 

OR INTI 001# I GOES #LATG# LONG 

19 CONTINUE 

C LONGITUDE# LATITUDE TO NIMBUS-ERB GRID 

CALL LLNIM{ XL ON# XL AT# NI MB# N IML AT# N IMLON ) 

PRINT1012 

PRINT1001#NIMB#NIMLAT#NIML0N 
C RETURN MIDPOINTS FROM GOES GRID 

IF(IERR.NE.O) G0T020 
CALL GOESLL(LATG#LONG#GLON#GLAT) 

PRINTIOIA 

PRINT1002#GL0N,GLAT 

20 CONTINUE 

C BACK TO LONGITUDE# LATITUDE FROM 3DNEPH GRID 

CALL NEPHLL(NEPHB# IH# NROW# NCOL# Z LON# ZLAT ) 

PRINT1004 
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PRINT1002>ZL0N> ZLAT 

C RETURN MIDPOINTS FROM ERBE GRIDS 
PRINT1016 
S-1.25 

DO 200 1-1,3 

$-$♦ 2 . 

CALL ERBLL(IROW(I)»ICOL(I)»S,CLON,CLAT) 

PPINT1003, S 
PRINT1009,CL0N,CLAT 
200 CONTINUE 

C RETURN MIDPOINTS FROM NIMBUS-ERB GRID 

CALL NIMLL(NIMB,CIMLON,CIMLAT) 

PRINT1013 

PRINT1002,CIMLON,CIMLAT 

GOTOl 

1000 FORMATC* 30NEPH LOCATION - BOX, ROW AND COL IDENTIFIERS*) 

1001 F0RMAT(3I6) 

1002 FORMATS* LONGITUDE AND LATITUDE ♦2F8.3) 

1003 FORMATdXFS.l,* DEGREE ERBE GRID*) 

lOOA F0RMAT(* return longitude, latitude FROM 3DNEPH BOX*) 

1005 F0RMAT(* northern HEMISPHERE*) 

1006 FORMATC* SOUTHERN HEMISPHERE*) 

1007 FORMAT(* LONG. OR LAT. OUTSIDE OF RANGE FOR GOES GRID*) 

1008 FORMAT(* GOES GRID BOX, ROW, COLUMN*) 

1009 FORMAT(* CENTER LONGITUDE AND LATITUDE *2F8.3) 

1010 FORMAT(//* INPUT LONGITUDE AND LATITUDE, FREE FORMAT*) 

1012 FORMAT!* NIMBUS-ERB TARGET AREA, LAT AND' LON REFERENCE*) 

1013 FORMAT!* RETURN LONGITUDE, LATITUDE FROM NIMBUS-ERB BOX*) 
101<> PORMAT!* RETURN LONGITUDE, LATITUDE FROM GOES GRID*) 

1015 FORMAT!* LOCATIONS ON 3 ERBE GRID SYSTEMS*) 

1016 FORMAT!* RETURN MIDPOINTS FOR 3 ERBE GRIDS*) 

1017 F0RMAT!1X2F10.3) 

9999 STOP 

END 
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SUBROUTINE LLNEPH<XLON>XLAT#NEPHB# 

CONVERTS LONGITUDE^ LATITUDE TO 3DNEPH 
BOX NUMBER AND ROW AND COLUMN INDICES. 

CONVERTS ANY LONGITUDE TO -180<-L0N6<»180 
ANGLE(X)=AM0D(X>360.) 

XL»ANGLE(XLON) 

IF{XL.GT.180.) XL-XL-360. 

IF(ABSdOO.-XL) .LT.l.E-^)XL-99.9 999 
IF(A3S(-80.-XL) .LT.l.E-1 )XL«-79.9999 

X»SIGN(1.,COSO( 350. +XL ) ) ♦SORT <6231 7. 6272* { 1 .-SIND ( ABS ( XL AT ) ) )/ 
C (l.+SIND( ABS(XLAT) > )/(l.+TAND(350.+XL)**2) ) 

Y»X*TAND(350.+XU 

I«257+INT(X+.5) 

J»257-INT(Y+.5) 

IFCXLAT.LT.O.) J-257+INT( Y+.5 > 

DO 100 K«li64 
Jl»(K-l)/8 
J1«J-J1*6A 
IF( J1.GT.64)GQT0100 
IR*K 
GOTOllO 
100 CONTINUE 
110 DO 120 LL-1>8 
IC-LL-1 
Il-I-IC+64 

IF(I1.GT.6^) GQT0120 
G0T0130 
120 CONTINUE 
130 NEPHB»IR+IC 
RETURN 
END 
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SUBROUTINE NEPHLL( IB#H> IRQW# ICOL, XLON, XLAT ) 
CONVERTS 3DNEPH BOX# HEMISPHERE# ROW# AND COLUMN 
BACK TO LONGITUDE AND LATITUDE 
INTEGER H#D 

ANGLE(X)«AM0D(X#360.)+180.-SIGN(180.»X) 

RO-57. 29577951308 
0-(IB-l)/8 

I=MOD( (IB-1)#8)*6^+IC0L 

J»D*64+IR0W 

X»I-257 

Y»257-J 

IF(H.E0.2) Y-J-257 
F-X^X+Y^Y 

F»( 62317. 6272-F)/ (623 17. 6272+F) 

IF (F.LT.DGOTOIO 
XLAT«90. 

G0T020 

XLAT«ASIN<F)*RO 
IF(H.EQ.2)XLAT»-XLAT 
TF(Y.NF.0.) G0T033 
IF(X.LT.0.)XL0N — 170. 

IF(X.EQ.O.) XL0N«0. 

IFCX.GT.O.) XLON-10. 

G0T034 

XL0N»ATAN2(Y#X)*RD 

XL0N»XL0N-350. 

XL0N=ANGLE(XL0N) 

RETURN 

END 


SUBROUTINE LLERB ( XLON# XL AT# S# IBOX# LAT# LON ) 

CONVERTS LONGITUDE# LATITUDE TO BOX AND ROW# COLUMN 
INDICATORS IN 2.5# 5# OR 10 DEGREE ERBE GRIDS. 
CONVERTS ANY LONGITUDE TO 0<-L0N6<360 
ANGLE (X)-AMOD(X# 360.) ♦180.-SIGN( 180.# X) 
N-3600/IFIX(S*10. ) 

XL«ANGLE(XLON) 

Nl-N/4+1 

LAT«MIN0{N/2#INT(FL0AT(N1)-XLAT/S+1.E-9) ) 
L0N-MIN0(N#INT(XL/S+l.*l.E-9) ) 

IBOX«(LAT-l) + N-*-LON 

RETURN 

END 


SUBROUTINE ERBLL ( L AT# LON# S# CLON# CL AT ) 

GIVES CENTER LATITUDE AND LONGITUDE FOR ANY ERB GRID 

CLAT»90.-S/2.-FL0AT(LAT-l)^S 

CLON» F LOAT ( LON-1 )*S+S/2. 

RETURN 

END 
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SUBROUTINE LLGOES ( XL0N» XLAT.> I GOES > LAT> LON/ I ERR ) 

C CONVERTS LONGITUDE TO BOXES IN A G0ES-LAN6LEY GRID. 

C CONVERTS ANY LONGITUDE TO -180<-L0NG<»180 . 

DIMENSION Z(3),D(3)>XLQ{3)#XHI(3) 

ANGLE (X)«AM0D<X^360. ) 

DATA Z/123.75#132.50f lAl.00/,0/2 .2 5^ 2.50> 3.00/ 

DATA XLn/-121.5>-130.0>-138.0/»XHI/-31.5>-30.0>-18.0/ 
I=INT(ABS(XLAT)/18. + 1. ) 

XL*ANGLE(XL0N) 

IF ( (XL.LT.XLOd ) ) .OR.(XL .GE.XHKI ) ) ) G0T020 

IF(ABS(XLAT).GT.45.)G0T020 

IERR-0 

LAT»MIN0(<t0, (INT( (47.25-XLAT)/2.25 + l.E-9) ) ) 

LON»MINO( AO>INT( (Z(I)+XL0N)/D(I)+l.E-9) ) 
IG0ES»(LAT-1)*A04-L0N 
GOTOlOO 
20 IERR»1 
100 RETURN 
END 


SUBROUTINE GOESLL ( L AT# LON. CLON,C L AT ) 

C CONVERTS GOES LAT#LON IDENTIFIERS TO CENTER LONG AND LAT. 
DIMENSION 0(3)»Z(3) 

DATA Z/123.75#132.50#1A1.00/# D/ 2 . 25# 2 . 50# 3 .00/ 

CL AT-^6.125-2.25*FL0AT(LAT) 

K»LAT 

IF(LAT.GT. 20) K-41-LAT 

I« INT ( FLOAT (28-K)/8.+l.E~9) 

CL0N=D(I)«FL0AT(L0N)-(Z(I)-0(I)/2.) 

RETURN 

end 
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SUBROUTINE LLNIMt XL0N> XL AT» IB>LAT#L0N) 

C CONVERTS LONGITUDE AND LATITUDE TO BOXES IN NIMBUS-ER8 GRID 

C CONVERTS ANY LONGITUDE TOO<«LON6<360 

DIMENSION IG(A0)»At20) 

COMMON/NIMDAT/IG# A 

DATA A/3.#9.>16.>20.>30.>36.#40.>'»5,>^8.j 3*60.>4*72.^ 4t80./ 
DATA IG/ 0» 3f 12, 28, 48> 78, 114> 154, 

199, 247, 307, 367, 427, 499, 571, 643, 

715, 795, 875, 955,1035,1115,1195,1275, 
1355,1427,1499,1571,1643,1703,1763,1823, 
1871,1916,1956, 1992,2022, 20«t2, 2058 , 2067 / 

ANGLE ( X) - AMOO (X, 360.) +180. -SIGN (1 80., X) 

XL»ANGLE(XLON) 

LAT«MIN0(40, (INT( (XLAT+90. ) /4. 5+1 . E-9 ) +1 ) ) 

I-LAT 

IF(LAT.GT.20) I-41-LAT 
C-360.-XL 
T«360./A<I) 

L0N»INT(C/T+l.E-9)+l 
IF(XL.EQ.O.) LON»LON-INT(A( I) ) 

IB«IG(LAT)+LON 

return 

END 


SUBROUTINE NI HLL ( I 8, C LON,C L AT ) 

C GIVES CENTER LONGITUDE AND LATITUDE FOR ANY NIMBUS-ERB BOX 

DIMENSION IG(40),A(20) 

COMMON/NIMDAT/IG, A 
00 10 I«l,40 
IF(IG(I) .LT.IB) GOTOlO 
LAT-I-1 
G0T020 
10 CONTINUE 
LAT»40 
20 J-LAT 

IF(J.GT.20) J-41-LAT 
LON-IB-IG(LAT) 

T«360./A( J) 

CLnN»360.-<T+FL0AT(L0N-l)+T/2.) 

CLAT — 90. +FL0AT(LAT-l)+4.5+2.2 5 
IF(CL0N.GT.180. )CL0N»CL0N-360. 

RETURN 

END 
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Sample Output From Program LLGRID 


INPUT LONGITUDE AND LATITUDE/ FREE FORMAT 

0.000 0.000 

3DNEPH LOCATION - BOX, ROW AND COL IDENTIFIERS 
NORTHERN HEMISPHERE 
^0 43 55 

LOCATIONS ON 3 ERBE GRID SYSTEMS 

2.5 DEGREE ERBE GRID 
5185 37 1 

5.0 DEGREE ERBE GRID 
1297 19 1 

10.0 DEGREE ERBE GRID 
325 10 1 

LONG. OR LAT. OUTSIDE OF RANGE FOR GOES GRID 
NIMBUS-ERB TARGET AREA, LAT AND LON REFERENCE 
1036 21 1 

RETURN LONGITUDE, LATITUDE FROM 3DNEPH BOX 
LONGITUDE AND LATITUDE .311 .017 

RETURN MIDPOINTS FOR 3 ERBE GRIDS 

2.5 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 1.250 -1.250 

5.0 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 2.500 -2.500 

10.0 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 5.000 -5.000 

RETURN LONGITUDE, LATITUDE FROM NIMBUS-ERB BOX 
LONGITUDE AND LATITUDE -2.250 2.250 


INPUT LONGITUDE AND LATITUDE, FREE FORMAT 
0.000 90.000 

3DNEPH LOCATION - BOX, ROW AND COL IDENTIFIERS 
NORTHERN HEMISPHERE 
37 1 1 

LOCATIONS ON 3 ERBE GRID SYSTEMS 

2.5 DEGREE ERBE GRID 
111 

5.0 DEGREE ERBE GRID 
111 

10.0 DEGREE ERBE GRID 
111 

LONG. OR LAT. OUTSIDE OF RANGE FOR GOES GRID 
NIMBUS-ERB TARGET AREA, LAT AND LON REFERENCE 
2068 40 1 

RETURN LONGITUDE, LATITUDE FROM 3DNEPH BOX 
LONGITUDE AND LATITUDE 0.000 90.000 

RETURN MIDPOINTS FOR 3 ERBE GRIDS 

2.5 DEGREE ERBE GRID 
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Sample Output From Program LLGRID - Continued 


CENTER LONGITUDE AND LATITUDE 

5.0 DEGREE ERBE GRID 
CENTER LONGITUDE AND LATITUDE 

10.0 DEGREE ERBE GRID 
CENTER LONGITUDE AND LATITUDE 
RETURN LONGITUDE, LATITUDE FROM 
LONGITUDE AND LATITUDE -60.000 


1.250 

08.750 

2.500 

87.500 

5.000 

85.000 

NIMBUS- 

ERB BOX 

87.750 



INPUT LONGITUDE AND LATITUDE, FREE FORMAT 
0.000 -90.000 

3DNEPH LOCATION - BOX, ROW AND COL IDENTIFIERS 
SOUTHERN HEMISPHERE 
37 1 1 

LOCATIONS ON 3 ERBE GRID SYSTEMS 
2.5 DEGREE ERBE GRID 
10225 72 1 

5.0 DEGREE ERBE GRID 
2521 36 1 

10.0 DEGREE ERBE GRID 
613 18 1 

LONG. OR LAT. OUTSIDE OF RANGE FOR GOES GRID 
NIMBUS-ERB TARGET AREA, LAT AND LON REFERENCE 
111 

RETURN LONGITUDE, LATITUDE FROM 3DNEPH BOX 
LONGITUDE AND LATITUDE 0.000 -90.000 
RETURN MIDPOINTS FOR 3 ERBE GRIDS 
2.5 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 1.250 -80.750 

5.0 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 2.500 -87.500 

10.0 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 5.000 -85.000 
RETURN LONGITUDE, LATITUDE FROM NIMBUS-ERB BOX 
LONGITUDE AND LATITUDE -60.000 -87.750 


INPUT LONGITUDE AND LATITUDE, FREE FORMAT 
-100.000 40.000 

3DNEPH LOCATION - BOX, ROW AND COL IDENTIFIERS 
NORTHERN HEMISPHERE 
44 45 26 

LOCATIONS ON 3 ERBE GRID SYSTEMS 


2.5 

DEGREE 

ERBE 

GRID 

2985 

21 

105 


5.0 

DEGREE 

ERBE 

GRID 

773 

11 

53 


10.0 

DEGREE 

ERBE 

GRID 

207 

6 

27 
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Sample Output From Program LLGRID - Concluded 


GOES GRID BOX, ROW# COLUMN 
93 3 13 

NIMBUS-ERB TARGET AREA, LAT AND LON REFERENCE 
1660 29 17 

RETURN LONGITUDE, LATITUDE FROM GOES GRID 
LONGITUDE AND LATITUDE -100.500 39.375 

RETURN LONGITUDE, LATITUDE FROM 3DNEPH BOX 
LONGITUDE AND LATITUDE 260.145 40.597 
RETURN MIDPOINTS FOR 3 ERBE GRIDS 

2.5 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 261.250 38.750 

5.0 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 262.500 37.500 

10.0 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 265.000 35.000 

RETURN LONGITUDE, LATITUDE FROM NiMBUS-ERB BOX 
LONGITUDE AND LATITUDE -99.000 38.250 


INPUT LONGITUDE AND LATITUDE, FREE FORMAT 
-75.000 -45.000 

3DNEPH LOCATION - BOX, ROW AND COL IDENTIFIERS 
SOUTHERN HEMISPHERE 
21 27 10 

LOCATIONS ON 3 ERBE GRID SYSTEMS 

2.5 DEGREE ERBE GRID 
7891 55 115 

5.0 DEGREE ERBE GRID 
2002 26 56 

10.0 DEGREE ERBE GRID 
497 14 29 

GOES GRID BOX, ROW, COLUMN 
1582 40 22 

NIMBUS-ERB TARGET AREA, LAT AND LON REFERENCE 
320 11 13 

RETURN LONGITUDE, LATITUDE FROM GOES GRID 
LONGITUDE AND LATITUDE -73.500 -43.875 
RETURN longitude, LATITUDE FROM 3DNEPH BOX 
LONGITUDE AND LATITUDE 285.042 -45.395 
RETURN MIDPOINTS FOR 3 ERBE GRIDS 

2.5 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 286.250 -46.250 

5.0 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 287.500 -47.500 

10.0 DEGREE ERBE GRID 

CENTER LONGITUDE AND LATITUDE 285.000 -45.000 
RETURN LONGITUDE, LATITUDE FROM NIMBUS-ERB BOX 
LONGITUDE AND LATITUDE -75.000 -42.750 
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PROGRAM LONGRA( INPUT# OUT PUT» TAPES -INPUT^ TAP E6-0UTPUT) 

C DETERMINES RIGHT ASCENSION FOR AN INPUT LONGITUDE AND 

C CALENDAR DATA. 

DIMENSION X(6) 

PRINTlOOl 
1 READ*# X,XLONG 

IF (X(1).EQ.99.) G0T0999 

PRINT100A#X#XL0NG 

CALL CALJULKX# WJD#FJD) 

XJD»«JD+FJD 

CALL SIDTIM(WJD»FJD»XLONG#RA) 

PRINT1002 

PRINT1003#X#XJD#XL0NG#RA 

GOTOl 

1001 FORMAT(//* INPUT CALENDAR DATE YEARtLAST 2 DIGITS )# MONTH# DAY# 
CHOUR#MINUTE#SEC# THEN L ONGI TUDE ( DEG ) . FREE FORMAT.*# 

C+ ENTER X(l)«99 TO STOP.*) 

1002 FDRMAT(* CALENDAR DATE# JULIAN DATE# LONGITUDE# RIGHT ASCENSION*) 

1003 F0RMAT(1X6F3.0,F15.6#2F8.2) 
lOOA FORMAT(1X6F4.0# F10.3) 

999 STOP 
END 


SUBROUTINE SIDTIM ( WJD# F J D# ALONG# ST ) 

CALCULATES LOCAL SIOERAL TIME 
A-JULIAN DATE AT 0 HRS. UNIVERSAL TIME 
B=NUMBEP OF MINUTES FROM 0 HRS. U.T. TO LOCAL TIME 
IF(FJD.LE..S) GO TO 1 
A-WJD+.S 
B-{FJD-.5)*14A0. 

60 TO 2 

1 A»WJD“.S 
B«(FJD+.5)*14A0. 

2 TU«(A-2415020.)/36525. 

TG0«99. 69098 33+36000. 7689*TU+. 0003 8708*TU**2+.2506B447*B+AL0NG 
ST«AM0D(TG0#360. ) +180. -SIGN (180.# TGO) 

RETURN 

END 
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SUBROUTINE C AL JUL 1 ( X, HJ D# F J 0 ) 

INPUT : X»YR(LAST TWO D I GITS ) , MONTH, D AY, HR, MIN, SEC> UT 
OUTPUT: WJD*WHOLE JULIAN DATE, F J D«FRACTI ONAL DATA 
DIMENSION X(6),A(12) 

D50-2A33282. 

YD»X( 1 1-A8. 

YL-YD/A. 

KYL»YL 

CK=KYL 

IF(YL-CK)1,1,3 
1 IF(X(2)-2.)4,4,3 
3 DS«CK 
G0T05 

A DS-CK-1. 

5 DS»DS+365.*( YD-2. ) 

DO 6 I«l,12 

6 A{I)»1. 

K«X(2) 

DO 7 I»K,12 

7 A(I)=0. 

DS«DS+31«*(A(1)+A(3)+A{5)+A(7)+A(8)+A(10)+A(12)) 

C+30.*(A{A)+A(6)+A(9)+A(11))+28.*A(2) 

DS»DS+X(3)-1, 

WND«DS 

FD»X( A)/2A.+X(5)/14A0.+X(6)/86A00. 

IF(F0-,A<39Qq99)9,8,e 

8 FJD-FO-,5 
WJO=l. 

GOTOlO 

9 FJD-FD+.5 
WJD-0. 

10 WJD-050+WJ0+WN0 
RETURN 
END 


Sample Output From Program LONGRA 


INPUT (UT) CALENDAR DATE YEAR(LAST 2 DIGITS ), MONTH, DAY, HOUR, MINUTE, SEC, 
THEN LONGITUDE(DEG) . FREE FORMAT, ENTER X(l)«99 TO STOP. 

81. 3. 21. 12, 0, 0. 0.000 

CALENDAR DATE, JULIAN DATE, LONGITUDE, RIGHT ASCENSION 
81. 3.21.12. 0. 0. 2AAA685. 000000 0.00 358.92 

81. 1. 1. 0. 0. 0. 0.000 

CALENDAR DATE, JULIAN DATE, LONGITUDE, RIGHT ASCENSION 
81. 1. 1. 0. 0, 0, 2AAA605. 500000 0.00 100,56 

81, 9. 29. 12, 0, 0. -75.000 

CALENDAR DATE, JULIAN DATE, LONGITUDE, RIGHT ASCENSION 
81. 9,29.12. 0. 0. 2AAA877. 000000 -75.00 113.16 

81. 9. 29. 13. 0. 0, -75,000 

CALENDAR DATE, JULIAN DATE, LONGITUDE, RIGHT ASCENSION 
81. 9.29.13. 0. 0, 2A4A877, 041667 -75.00 128.20 
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PROGRAM ZENITH<INPUT#0UTPUT^TAPE5»INPUT>TAPE6«0UTPUT) 

C CALCULATES COSINE OF SOLAR ZENITH ANGLE FOR INPUT OF 

C LONGITUDE, LATITUDE, GMT CALENDAR DATE. 

DIMENSION X(6) 

PRINTlOOl 

1 READ+,X,XLON,XLAT 

IF(X(l).EQ.O.) G0TQ999 
PRINT100A,X,XL0N,XLAT 
CALL CALJUL1(X,WJD»FJD) 

XJD»WJD+F JD 

CALL r.ZEN(XJD,XLAT,XLON»SOLZEN) 

CSZEN«COSD(SOLZEN) 

PRINT1002 

PRINTlO03fX,XJD,XLON,XLAT,SOLZEN,CSZEN 

GOTOl 

1001 FORMAT(//* INPUT YR(LAST 2 DIGITS), MON, DAY, HR, MIN, SEC , LONG, LAT*/ 
C* SET X(1)«0 TO TERMINATE^) 

1002 FORMAT(* CAL., JUL. DATE, LONG, LAT, SOLAR ZENITH ANGLE £ COSt) 

1003 F0RMAT(lX6F3.0,F15.5,2F7.2,F10.4,F7.<r) 
lOOA F0RMAT(1X6FA.0,2F10.3) 

999 STOP 
END 


subroutine CZEN(XJ0,XLAT,XL0N,S0LZEN) 

C finds COSINE OF SOLAR ZENITH ANGLE FOR JULIAN DATE, 

C LATITUDE, AND LONGITUDE 

CALL EEARTH(XJD,X,Y,Z,OX,DY,DZ) 

P=SQRT(X*X+Y*YfZ«Z) 

XS»-X/R 
YS— Y/R 
ZS—Z/R 

CALL RECEQ(XJD,XS,YS,ZS,XQ, YQ,ZQ) 

WJD«INT(XJD) 

FJD-XJD-WJD 

CALL SIDTIM(WJD,FJD,XLON,RA) 

XE»COSD(RA)*COSD(XLAT) 

YE-SIND(RA)*C0SD(XLAT) 

ZE«SIND(XLAT) 

C COSZEN is dot PRODUCT OF UNIT SUN AND POSITION VECTORS 
SOLZEN-ACOS (XE+XQ+ YE* YQ+ZE+ZQ )/. 01745329252 
RETURN 
END 


SUBROUTINE C ALJ UL 1 ( X, W J D , F J 0) 

(See 6ubn.0iUUnz LiAtingi ^on. pAogAom LOMGKA) 
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SUBROUTINE EEARTH( JD# XHE#YHE> ZHE> DXHE> DYHE^DZHE) 

THIS SUBROUTINE COMPUTES THE HELIOCENTRIC POSITION AND VELOCITY OF 
THE EARTH IN MEAN EQUINOX AND ECLIPTIC OF DATE COORDINATE SYSTEM. 
THIS ROUTINE CALLS SUBROUTINES TINVS AND CONCAR, 

JD - JULIAN DATE 
XHE>YHE>ZHE - POSITION OF EARTH 
DXHE>DYHE>DZHE - VELOCITY OF EARTH 

REAL JD 

ANGLE ( X) « AM0D( X> 360. ) +18 0.-SlGN(180.,X) 

AU=«1495988A5. S USUN*1 . 3271 AilE + 11 
DR=.017A532925199A3 
R0»57. 2957795130823 
C 

D»J0-2415020. 

CD-D/10000. 

TE-D/36525. 

C 

AE-1.00000023+AU 

EE*0.0167510A-O.OOOOA180+TE~0.000000126*TE+*2 

XIE-0.0 

WE»101.220833+0.0000A7068*D+0.0000339*C0^+2+0.00000007*CD+^3 

0E«0.0 

XME«ANGLE{358.A75845+0.985600267+D-0. 00001 12*CD**2-0.00000007*CD»* 
13) 

C 

CALL TINVS( XME*DR^ EE,ECE#FE ) 

CALL CONCAR(AE#EE#XIE>WE>OE>FE*RD#XHE»YHEiZHE,DXHE»DYHE>DZHE»USUN) 
C 

RETURN 

END 
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SUBROUTINE C ONC AR ( A, E» X I # W> 0> F# X» Y# Z» DX> DY> DZ#U ) 

note: this is A SINGLE PRECISION VERSION OF CONCAR 

0R«.01745329251994f3 

FR aOR*F 

WFR«OR*(W+F) 

0R»DR*0 

XIR»0R*XI 

R»A*(l.-E*E)/(l.+E*COS(FR)) 

V»SQRT(U*(2./R-1./A) ) 

GAM-ATAN{E*SIN(FR)/a. + EKOS{FR) ) ) 

WFGR=*WFR~GAH 

CUF«COS(WFR) 

SUF-SIN(WFR) 

SQ»SINtOR) 

CO»COS(OR) 

SI-SIN(XIR) 

CI-COS(XIR) 

SWFG»SIN(WFGR) 

CWFG«COS(WFGR) 

X»R^(CWF*CO-SWF^SO*CI) 

Y»R*(CWF*SO+SWF*COtCI) 

Z«R*(SWF + SI ) 

DX»V* (-SWFG*CO-CWFG«SO*CI) 
DY-V+(-SWFG*SO+CWFGtCO*CI) 

DZ-V^(CWFG^SI) 

RETURN 

END 


SUBROUTINE TINVS ( M, E> EC >F ) 

REAL M,MO 

DATA PI/3.141592653589793/ 

ASINH(X) -SIGN (ALOG(ABS( X)+ SORT (X**2+l. ) ),X) 

IFIE.GE.l. )GO TO 100 

EC«M 

10 MO«EC-E*SIN(EC) 

DM=M-MO 

DE »DM/ (l.-E*COS{EC) ) 

EC-EC+DE 

IF(ABS(DE) .GT.l.E-12 ) GO TO 10 

HEC» EC/2. 

HF-ATAN(SQRT( (1.+E)/(1.~E) )*SIN(HEC)/C0S(HEC) ) 
IFCHF.LT.O. )HF»HF<-PI 
F«2.*HF 
GO TO 800 

100 CONTINUE 
6C»ASINH(M/E ) 

101 MO«E+SINH(EC)-EC 
DM=M-MO 

DE-DN/(E*C0SH(EC)-1.) 

EC»EC+DE 

IF(A3S(DE). GT.l.E-12 ) GO TO 101 

F»2. tATANI SORT (( E+1.0 )/( E-1.0) )*TANH{ EC/2.0 ) ) 
800 RETURN 
END 



oooooooooo 


APPENDIX 


SUBROUTINE REC EQ ( JD > X EC # YEC> ZEC» XEQf YEQ»ZEQ) 

THIS SUBROUTINE ROTATES A VECTOR FROM GEOCENTRIC# ECLIPTIC# TO 
THE GEOCENTRIC# EARTH EQUATORIAL COORDINATE SYSTEM 

JO - JULIAN DATE 

XEC#YEC#ZEC - COMPONENTS OF THE VECTOR IN THE GEOCENTRIC# ECLIPTIC 
COORDINATE SYSTEM 

XEQ#YEQ#ZEQ - COMPONENTS OF THE VECTOR IN THE GEOCENTRIC# EARTH 
EQUATORIAL# COORDINATE SYSTEM 

REAL JD 

DR-.017A53292519943 
TE»( JD-241 5020,) 736525 . 

XIE«23.452294-0.0130125*TE-0.00000164*TE*+2+0.000000503*TE**3 

C»COS(XIE+DR) 

S*SIN(XIE*DR) 

XEQ-XEC 

YEQ«YEC*C-ZEC+S 

ZEQ*YEC*S+ZEC*C 

RETURN 

END 


SUBROUTINE SIDTIMC WJD# FJD# ALONG# ST) 

(See ^ub/LoivUne. tlitingi {^on. ph-ognam LONGRA} 


Sample Output From Program ZENITH 


INPUT YR(LAST 2 DIGITS)# MON# DAY#HR#MIN#SEC# LONG#LAT 
SET X (1 )=0 TO TERMINATE 
81. 9. 29. 12. 0. 0. 0.000 37.000 

CAL.# JUL. DATE# LONG# LAT# SOLAR ZENITH ANGLE & COS 
81. 9.29,12. 0. 0, 2444877.00000 0.00 37.00 39.5511 .7711 

81. 9. 29. 13. 0. 0. -75.000 37,000 

CAL.# JUL. DATE# LONG# LAT# SOLAR ZENITH ANGLE C COS 
81. 9.29.13. 0. 0. 2444877.04167 -75.00 37.00 66.3318 ,4014 

81. 3. 21. 12. 0. 0. 0.000 0,000 

CAL.# JUL. DATE# LONG# LAT# SOLAR ZENITH ANGLE £ COS 
81. 3.21.12. 0. 0. 2444685.00000 0.00 0.00 1.8376 .9995 
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Figure 2.- The Langley GOES grid system. 
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Figure 5.- Relationship between longitude and right ascension 
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Figure 6.- Definition of solar right ascension, 
declination, and zenith angles in terms of 
solar and ground observer unit position vector. 
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